home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_temme.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  6.4 KB  |  220 lines
  1. /* specfunc/bessel_temme.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. /* Calculate series for Y_nu and K_nu for small x and nu.
  23.  * This is applicable for x < 2 and |nu|<=1/2.
  24.  * These functions assume x > 0.
  25.  */
  26. #include <config.h>
  27. #include <gsl/gsl_math.h>
  28. #include <gsl/gsl_errno.h>
  29. #include <gsl/gsl_mode.h>
  30. #include "bessel_temme.h"
  31.  
  32. #include "chebyshev.h"
  33. #include "cheb_eval.c"
  34.  
  35. /* nu = (x+1)/4, -1<x<1, 1/(2nu)(1/Gamma[1-nu]-1/Gamma[1+nu]) */
  36. static double g1_dat[14] = {
  37.   -1.14516408366268311786898152867,
  38.    0.00636085311347084238122955495,
  39.    0.00186245193007206848934643657,
  40.    0.000152833085873453507081227824,
  41.    0.000017017464011802038795324732,
  42.   -6.4597502923347254354668326451e-07,
  43.   -5.1819848432519380894104312968e-08,
  44.    4.5189092894858183051123180797e-10,
  45.    3.2433227371020873043666259180e-11,
  46.    6.8309434024947522875432400828e-13,
  47.    2.8353502755172101513119628130e-14,
  48.   -7.9883905769323592875638087541e-16,
  49.   -3.3726677300771949833341213457e-17,
  50.   -3.6586334809210520744054437104e-20
  51. };
  52. static cheb_series g1_cs = {
  53.   g1_dat,
  54.   13,
  55.   -1, 1,
  56.   7
  57. };
  58.  
  59. /* nu = (x+1)/4, -1<x<1,  1/2 (1/Gamma[1-nu]+1/Gamma[1+nu]) */
  60. static double g2_dat[15] = 
  61. {
  62.   1.882645524949671835019616975350,
  63.  -0.077490658396167518329547945212,  
  64.  -0.018256714847324929419579340950,
  65.   0.0006338030209074895795923971731,
  66.   0.0000762290543508729021194461175,
  67.  -9.5501647561720443519853993526e-07,
  68.  -8.8927268107886351912431512955e-08,
  69.  -1.9521334772319613740511880132e-09,
  70.  -9.4003052735885162111769579771e-11,
  71.   4.6875133849532393179290879101e-12,
  72.   2.2658535746925759582447545145e-13,
  73.  -1.1725509698488015111878735251e-15,
  74.  -7.0441338200245222530843155877e-17,
  75.  -2.4377878310107693650659740228e-18,
  76.  -7.5225243218253901727164675011e-20
  77. };
  78. static cheb_series g2_cs = {
  79.   g2_dat,
  80.   14,
  81.   -1, 1,
  82.   8
  83. };
  84.  
  85.  
  86. static
  87. int
  88. gsl_sf_temme_gamma(const double nu, double * g_1pnu, double * g_1mnu, double * g1, double * g2)
  89. {
  90.   const double anu = fabs(nu);    /* functions are even */
  91.   const double x = 4.0*anu - 1.0;
  92.   gsl_sf_result r_g1;
  93.   gsl_sf_result r_g2;
  94.   cheb_eval_e(&g1_cs, x, &r_g1);
  95.   cheb_eval_e(&g2_cs, x, &r_g2);
  96.   *g1 = r_g1.val;
  97.   *g2 = r_g2.val;
  98.   *g_1mnu = 1.0/(r_g2.val + nu * r_g1.val);
  99.   *g_1pnu = 1.0/(r_g2.val - nu * r_g1.val);
  100.   return GSL_SUCCESS;
  101. }
  102.  
  103.  
  104. int
  105. gsl_sf_bessel_Y_temme(const double nu, const double x,
  106.                       gsl_sf_result * Ynu,
  107.                       gsl_sf_result * Ynup1)
  108. {
  109.   const int max_iter = 15000;
  110.   
  111.   const double half_x = 0.5 * x;
  112.   const double ln_half_x = log(half_x);
  113.   const double half_x_nu = exp(nu*ln_half_x);
  114.   const double pi_nu   = M_PI * nu;
  115.   const double alpha   = pi_nu / 2.0;
  116.   const double sigma   = -nu * ln_half_x;
  117.   const double sinrat  = (fabs(pi_nu) < GSL_DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));
  118.   const double sinhrat = (fabs(sigma) < GSL_DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);
  119.   const double sinhalf = (fabs(alpha) < GSL_DBL_EPSILON ? 1.0 : sin(alpha)/alpha);
  120.   const double sin_sqr = nu*M_PI*M_PI*0.5 * sinhalf*sinhalf;
  121.   
  122.   double sum0, sum1;
  123.   double fk, pk, qk, hk, ck;
  124.   int k = 0;
  125.   int stat_iter;
  126.  
  127.   double g_1pnu, g_1mnu, g1, g2;
  128.   int stat_g = gsl_sf_temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);
  129.  
  130.   fk = 2.0/M_PI * sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);
  131.   pk = 1.0/M_PI /half_x_nu * g_1pnu;
  132.   qk = 1.0/M_PI *half_x_nu * g_1mnu;
  133.   hk = pk;
  134.   ck = 1.0;
  135.  
  136.   sum0 = fk + sin_sqr * qk;
  137.   sum1 = pk;
  138.  
  139.   while(k < max_iter) {
  140.     double del0;
  141.     double del1;
  142.     double gk;
  143.     k++;
  144.     fk  = (k*fk + pk + qk)/(k*k-nu*nu);
  145.     ck *= -half_x*half_x/k;
  146.     pk /= (k - nu);
  147.     qk /= (k + nu);
  148.     gk  = fk + sin_sqr * qk;
  149.     hk  = -k*gk + pk; 
  150.     del0 = ck * gk;
  151.     del1 = ck * hk;
  152.     sum0 += del0;
  153.     sum1 += del1;
  154.     if(fabs(del0) < 0.5*(1.0 + fabs(sum0))*GSL_DBL_EPSILON) break;
  155.   }
  156.  
  157.   Ynu->val   = -sum0;
  158.   Ynu->err   = (2.0 + 0.5*k) * GSL_DBL_EPSILON * fabs(Ynu->val);
  159.   Ynup1->val = -sum1 * 2.0/x;
  160.   Ynup1->err = (2.0 + 0.5*k) * GSL_DBL_EPSILON * fabs(Ynup1->val);
  161.  
  162.   stat_iter = ( k >= max_iter ? GSL_EMAXITER : GSL_SUCCESS );
  163.   return GSL_ERROR_SELECT_2(stat_iter, stat_g);
  164. }
  165.  
  166.  
  167. int
  168. gsl_sf_bessel_K_scaled_temme(const double nu, const double x,
  169.                              double * K_nu, double * K_nup1, double * Kp_nu)
  170. {
  171.   const int max_iter = 15000;
  172.  
  173.   const double half_x    = 0.5 * x;
  174.   const double ln_half_x = log(half_x);
  175.   const double half_x_nu = exp(nu*ln_half_x);
  176.   const double pi_nu   = M_PI * nu;
  177.   const double sigma   = -nu * ln_half_x;
  178.   const double sinrat  = (fabs(pi_nu) < GSL_DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));
  179.   const double sinhrat = (fabs(sigma) < GSL_DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);
  180.   const double ex = exp(x);
  181.  
  182.   double sum0, sum1;
  183.   double fk, pk, qk, hk, ck;
  184.   int k = 0;
  185.   int stat_iter;
  186.  
  187.   double g_1pnu, g_1mnu, g1, g2;
  188.   int stat_g = gsl_sf_temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);
  189.  
  190.   fk = sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);
  191.   pk = 0.5/half_x_nu * g_1pnu;
  192.   qk = 0.5*half_x_nu * g_1mnu;
  193.   hk = pk;
  194.   ck = 1.0;
  195.   sum0 = fk;
  196.   sum1 = hk;
  197.   while(k < max_iter) {
  198.     double del0;
  199.     double del1;
  200.     k++;
  201.     fk  = (k*fk + pk + qk)/(k*k-nu*nu);
  202.     ck *= half_x*half_x/k;
  203.     pk /= (k - nu);
  204.     qk /= (k + nu);
  205.     hk  = -k*fk + pk;
  206.     del0 = ck * fk;
  207.     del1 = ck * hk;
  208.     sum0 += del0;
  209.     sum1 += del1;
  210.     if(fabs(del0) < 0.5*fabs(sum0)*GSL_DBL_EPSILON) break;
  211.   }
  212.   
  213.   *K_nu   = sum0 * ex;
  214.   *K_nup1 = sum1 * 2.0/x * ex;
  215.   *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
  216.  
  217.   stat_iter = ( k == max_iter ? GSL_EMAXITER : GSL_SUCCESS );
  218.   return GSL_ERROR_SELECT_2(stat_iter, stat_g);
  219. }
  220.